/* adt_complex.c
 * A complex (single precision only), bicomplex and quaternion numbers and arithmetic.
 * ADT (Altai Diffraction Team)
 *
 * Copyright (C) 2010 Vsevolod Scherbinin
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */
#include <math.h>
// #include "adt_complex.h"

#if 0
#if (ONHOST)
#define ADTFUNTYPE 
#define ADTFUNNAME(fn) adt_fn
#endif
#if (ONDEVICE)
#define ADTFUNTYPE  __device__
#define ADTFUNNAME(fn) adt_cuda_fn
#endif
#endif

// #if (ONDEVICE)
// #define ADTFUNTYPE  __device__
// #define ADTFUNNAME(fn) adt_cuda_complex_##fn
// #define ADTFUNNAMEB(fn) adt_cuda_bicomplex_##fn
// #else
// #define ADTFUNTYPE 
// #define ADTFUNNAME(fn) adt_complex_##fn
// #define ADTFUNNAMEB(fn) adt_bicomplex_##fn
// #endif
// #define ADTPRECN(fn) fn##_float


//returning complex number x+iy
ADTFUNTYPE adt_complex_float
ADTFUNNAME(rect) (float x, float y)
{
  adt_complex_float w;
  w.dat[0]=x;
  w.dat[1]=y;
  return w;
}

//returning complex number r \cos theta + i r \sin theta
ADTFUNTYPE adt_complex_float
ADTFUNNAME(polar) (float r, float theta)
{
  adt_complex_float w;
  w.dat[0] = r*cosf(theta);
  w.dat[1] = r*sinf(theta);
  return w;
}

//This function returns the argument of the complex number z, \arg(z), where -\pi < \arg(z) <= \pi. 
ADTFUNTYPE float
ADTFUNNAME(arg) (adt_complex_float z)
{
  float w;
  w = atan2f (z.dat[1], z.dat[0]);
  return w;
}

//This function returns the magnitude of the complex number z, |z|. 
ADTFUNTYPE float
ADTFUNNAME(abs) (adt_complex_float z)
{
  float w;
  w = sqrtf(z.dat[0]*z.dat[0] + z.dat[1]*z.dat[1]);
  return w;
}

//This function returns the squared magnitude of the complex number z, |z|^2. 
ADTFUNTYPE float
ADTFUNNAME(abs2) (adt_complex_float z)
{
  float w;
  w = z.dat[0]*z.dat[0] + z.dat[1]*z.dat[1];
  return w;
}

//This function returns the sum of the complex numbers z1 and z2, w=z1+z2. 
ADTFUNTYPE adt_complex_float
ADTFUNNAME(add) (adt_complex_float z1, adt_complex_float z2)
{
  adt_complex_float w;
  w.dat[0] = z1.dat[0] + z2.dat[0];
  w.dat[1] = z1.dat[1] + z2.dat[1]; 
  return w;
}

//This function returns the difference of the complex numbers z1 and z2, w=z1-z2. 
ADTFUNTYPE adt_complex_float
ADTFUNNAME(sub) (adt_complex_float z1, adt_complex_float z2)
{
  adt_complex_float w;
  w.dat[0] = z1.dat[0] - z2.dat[0];
  w.dat[1] = z1.dat[1] - z2.dat[1]; 
  return w;
}

//This function returns the product of the complex numbers z1 and z2, w=z1z2. 
ADTFUNTYPE adt_complex_float
ADTFUNNAME(mul) (adt_complex_float z1, adt_complex_float z2)
{
  adt_complex_float w;
  w.dat[0] = z1.dat[0] * z2.dat[0] - z1.dat[1] * z2.dat[1];
  w.dat[1] = z1.dat[0] * z2.dat[1] + z1.dat[1] * z2.dat[0]; 
  return w;
}

//This function returns the quotient of the complex numbers z1 and z2, w=z1/z2. 
ADTFUNTYPE adt_complex_float
ADTFUNNAME(div) (adt_complex_float z1, adt_complex_float z2)
{
  float den_div;
  adt_complex_float w;
  den_div = z2.dat[0]*z2.dat[0] + z2.dat[1]*z2.dat[1];
  w.dat[0] = z1.dat[0] * z2.dat[0] + z1.dat[1] * z2.dat[1];
  w.dat[1] = z1.dat[1] * z2.dat[0] - z1.dat[0] * z2.dat[1];
  w.dat[0] = w.dat[0]/den_div;
  w.dat[1] = w.dat[1]/den_div;
  return w;
}

//This function returns the sum of the complex number z and the real number x, w=z+x. 
ADTFUNTYPE adt_complex_float
ADTFUNNAME(add_real) (adt_complex_float z, float x)
{
  adt_complex_float w;
  w.dat[0] = z.dat[0] + x;
  w.dat[1] = z.dat[1]; 
  return w;
}

//This function returns the difference of the complex number z and the real number x, w=z-x. 
ADTFUNTYPE adt_complex_float
ADTFUNNAME(sub_real) (adt_complex_float z, float x)
{
  adt_complex_float w;
  w.dat[0] = z.dat[0] - x;
  w.dat[1] = z.dat[1]; 
  return w;
}

//This function returns the product of the complex number z and the real number x, w=zx. 
ADTFUNTYPE adt_complex_float
ADTFUNNAME(mul_real) (adt_complex_float z, float x)
{
  adt_complex_float w;
  w.dat[0] = z.dat[0] * x;
  w.dat[1] = z.dat[1] * x; 
  return w;
}

//This function returns the quotient of the complex number z and the real number x, w=z/x. 
ADTFUNTYPE adt_complex_float
ADTFUNNAME(div_real) (adt_complex_float z, float x)
{ 
  adt_complex_float w;
  w.dat[0] = z.dat[0]/x;
  w.dat[1] = z.dat[1]/x; 
  return w;
}

//This function returns the sum of the complex number z and the imaginary number iy, w=z+iy. 
ADTFUNTYPE adt_complex_float
ADTFUNNAME(add_imag) (adt_complex_float z, float y)
{
  adt_complex_float w;
  w.dat[0] = z.dat[0];
  w.dat[1] = z.dat[1] + y; 
  return w;
}

//This function returns the difference of the complex number z and the imaginary number iy, w=z-iy. 
ADTFUNTYPE adt_complex_float
ADTFUNNAME(sub_imag) (adt_complex_float z, float y)
{
  adt_complex_float w;
  w.dat[0] = z.dat[0];
  w.dat[1] = z.dat[1] - y; 
  return w;
}

//This function returns the product of the complex number z and the imaginary number iy, w=z*(iy).
ADTFUNTYPE adt_complex_float
ADTFUNNAME(mul_imag) (adt_complex_float z, float y)
{
  adt_complex_float w;
  w.dat[1] = z.dat[0] * y;
  w.dat[0] = - z.dat[1] * y; 
  return w;
}

//This function returns the quotient of the complex number z and the imaginary number iy, w=z/(iy).
ADTFUNTYPE adt_complex_float
ADTFUNNAME(div_imag) (adt_complex_float z, float y)
{ 
  adt_complex_float w;
  w.dat[0] = z.dat[1]/y;
  w.dat[1] = -z.dat[0]/y; 
  return w;
}

//This function returns the complex conjugate of the complex number z, z^* = x - i y. 
ADTFUNTYPE adt_complex_float
ADTFUNNAME(conjugate) (adt_complex_float z)
{
  adt_complex_float w;
  w.dat[0] = z.dat[0];
  w.dat[1] = -z.dat[1];
  return w;
}

//This function returns the inverse, or reciprocal, of the complex number z, 1/z = (x - i y)/(x^2 + y^2). 
ADTFUNTYPE adt_complex_float
ADTFUNNAME(inverse) (adt_complex_float z)
{
  adt_complex_float w;
  float den;
  den = z.dat[0]*z.dat[0] + z.dat[1]*z.dat[1];
  w.dat[0] = z.dat[0]/den;
  w.dat[1] = -z.dat[1]/den;
  return w;
}

//This function returns the negative of the complex number z, -z = (-x) + i(-y). 
ADTFUNTYPE adt_complex_float
ADTFUNNAME(negative) (adt_complex_float z)
{
  adt_complex_float w;
  w.dat[0] = -z.dat[0];
  w.dat[1] = -z.dat[1];
  return w;
}

//This function returns the square root of the complex number z, \sqrt z. The branch cut is the negative real axis. The result always lies in the right half of the complex plane.
ADTFUNTYPE adt_complex_float
ADTFUNNAME(sqrt) (adt_complex_float z)
{
  adt_complex_float w;
  float csqrt_arg, csqrt_abs;
  csqrt_abs = sqrtf(ADTFUNNAME(abs)(z));
  csqrt_arg = .5*ADTFUNNAME(arg)(z);
  w.dat[0] = csqrt_abs*cosf(csqrt_arg);
  w.dat[1] = csqrt_abs*sinf(csqrt_arg);
  return w;
}

//This function returns the complex exponential of the complex number z, \exp(z).
ADTFUNTYPE adt_complex_float
ADTFUNNAME(exp) (adt_complex_float z)
{
  adt_complex_float w;
  float reexp;
  reexp = exp(z.dat[0]);
  w.dat[0] = reexp*cosf(z.dat[1]);
  w.dat[1] = reexp*sinf(z.dat[1]);
  return w;
}



//Unary function for bicomplex type
//returning bicomplex number x+jy+iu+ijv
ADTFUNTYPE adt_bicomplex_float
ADTFUNNAMEB(rect) (float x, float y, float u, float v)
{
  adt_bicomplex_float b;
  b.dat[0] = x;
  b.dat[1] = y;
  b.dat[2] = u;
  b.dat[3] = v;
  return b;
}

//returning bicomplex number Re(z1)+jIm(z1)+iRe(z2)+ijIm(z2)
ADTFUNTYPE adt_bicomplex_float
ADTFUNNAMEB(cmplx) (adt_complex_float z1, adt_complex_float z2)
{
  adt_bicomplex_float b;
  b.dat[0] = z1.dat[0];
  b.dat[1] = z1.dat[1];
  b.dat[2] = z2.dat[0];
  b.dat[3] = z2.dat[1];
  return b;
}

//returning the two first components as bicomplex number
ADTFUNTYPE adt_complex_float
ADTFUNNAMEB(real) (adt_bicomplex_float b)
{
  adt_complex_float z;
  z.dat[0] = b.dat[0];
  z.dat[1] = b.dat[1];
  return z;
}

//returning the two last components as bicomplex number
ADTFUNTYPE adt_complex_float
ADTFUNNAMEB(imag) (adt_bicomplex_float b)
{
  adt_complex_float z;
  z.dat[0] = b.dat[2];
  z.dat[1] = b.dat[3];
  return z;
}

//This function returns the magnitude of the bicomplex number b, |b|.
ADTFUNTYPE float
ADTFUNNAMEB(abs) (adt_bicomplex_float z)
{
  float a1, a2;
  a1 = z.dat[0]*z.dat[0] - z.dat[1]*z.dat[1] + z.dat[2]*z.dat[2] - z.dat[3]*z.dat[3];
  a1 = a1*a1;
  a2 = z.dat[0]*z.dat[1] + z.dat[2]*z.dat[3];
  a2 = 4.*a2*a2;
  return sqrtf(sqrtf(a1+a2));
}

//This function returns the quadrupled magnitude of the bicomplex number b, |b|^4.
ADTFUNTYPE float
ADTFUNNAMEB(abs4) (adt_bicomplex_float z)
{
  float a1, a2;
  a1 = z.dat[0]*z.dat[0] - z.dat[1]*z.dat[1] + z.dat[2]*z.dat[2] - z.dat[3]*z.dat[3];
  a1 = a1*a1;
  a2 = z.dat[0]*z.dat[1] + z.dat[2]*z.dat[3];
  a2 = 4.*a2*a2;
  return a1+a2;
}

//This function returns the sum of the complex numbers b1 and b2, b=b1+b2. 
ADTFUNTYPE adt_bicomplex_float
ADTFUNNAMEB(add) (adt_bicomplex_float b1, adt_bicomplex_float b2)
{
  adt_bicomplex_float w;
  w.dat[0] = b1.dat[0] + b2.dat[0];
  w.dat[1] = b1.dat[1] + b2.dat[1]; 
  w.dat[2] = b1.dat[2] + b2.dat[2]; 
  w.dat[3] = b1.dat[3] + b2.dat[3]; 
  return w;
}

//This function returns the difference of the bicomplex numbers b1 and b2, b=b1-b2. 
ADTFUNTYPE adt_bicomplex_float
ADTFUNNAMEB(sub) (adt_bicomplex_float b1, adt_bicomplex_float b2)
{
  adt_bicomplex_float w;
  w.dat[0] = b1.dat[0] - b2.dat[0];
  w.dat[1] = b1.dat[1] - b2.dat[1]; 
  w.dat[2] = b1.dat[2] - b2.dat[2]; 
  w.dat[3] = b1.dat[3] - b2.dat[3]; 
  return w;
}

//This function returns the product of the bicomplex numbers b1 and b2, b=b1b2. 
/* a b = (a_0 b_0 - a_1 b_1 - a_2 b_2 + a_3 b_3) + j (a_0 b_1 + a_1 b_0 - a_2 b_3 - a_3 b_2) + 
+ i (a_2 b_0 - a_3 b_1 - a_0 b_2 + a_2 b_3) + ij (a_2 b_1 + a_3 b_0 + a_0 b_3 + a_2 b_2)*/


ADTFUNTYPE adt_bicomplex_float
ADTFUNNAMEB(mul) (adt_bicomplex_float b1, adt_bicomplex_float b2)
{
  adt_bicomplex_float w;
  w.dat[0] = b1.dat[0]*b2.dat[0] - b1.dat[1]*b2.dat[1] - b1.dat[2]*b2.dat[2] - b1.dat[3]*b2.dat[3];
  w.dat[1] = b1.dat[1]*b2.dat[0] + b1.dat[0]*b2.dat[1] - b1.dat[3]*b2.dat[2] - b1.dat[2]*b2.dat[3];
  w.dat[2] = b1.dat[2]*b2.dat[0] - b1.dat[3]*b2.dat[1] - b1.dat[1]*b2.dat[3] + b1.dat[0]*b2.dat[2];
  w.dat[3] = b1.dat[0]*b2.dat[3] + b1.dat[1]*b2.dat[2] - b1.dat[2]*b2.dat[1] - b1.dat[3]*b2.dat[0];
  return w;
}

//This function returns the product of the bicomplex number b and the real number x, w=bx.
ADTFUNTYPE adt_bicomplex_float
ADTFUNNAMEB(mul_real) (adt_bicomplex_float b, float x)
{
  adt_bicomplex_float w;
  w.dat[0] = b.dat[0]*x;
  w.dat[1] = b.dat[1]*x; 
  w.dat[2] = b.dat[2]*x; 
  w.dat[3] = b.dat[3]*x; 
  return w;
}

//This function returns the quotient of the bicomplex number b and the real number x, w=b/x.
ADTFUNTYPE adt_bicomplex_float
ADTFUNNAMEB(div_real) (adt_bicomplex_float b, float x)
{
  adt_bicomplex_float w;
  w.dat[0] = b.dat[0]/x;
  w.dat[1] = b.dat[1]/x; 
  w.dat[2] = b.dat[2]/x; 
  w.dat[3] = b.dat[3]/x; 
  return w;
}

#undef ADTFUNTYPE 
#undef ADTFUNNAME